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In  this  report  we  list  the  Fortran  subroutines  for 
evaluating  the  confluent  hypergeometric  functions  M(a,b:x)  and 
U(a,b;x).  These  subroutines  use  the  stable  recurrence 
relations  given  e.g.  in  Wimp. 
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I ntroduct i on 


It  is  well  known  that  the  ordinary  differential  equation 


d2v 


dy 


x  =-%  +  (1-x)  -  ay  =  0 

dx2  dx  J 


has  a  solution 


y(x)  =  AM(a, 1 ;x)  +  BU(a,l;x) 


if  a  is  not  a  negative  integer. 

This  problem  arises  e.g.  when  solving  the  linearized 
shallow  water  equations  with  the  full  linear  variation  in  depth 
included  (see  Williams,  Stan iforth  and  Neta,  [1] ) . 

The  computation  of  the  confluent  hypergeometric  functions 
is  based  on  the  Miller  algorithm  (see  e.g.  Wimp.  [2]).  In 
general ,  one  has  a  second  order  difference  equation 

z(n)  +  a ( n ) z ( n  +  1 )  +  b(n)z(n+2)  =  0.  n  >  0,  b(n)  ^  0  . 

If  b(n)  =  0  for  some  n,  in  some  cases  one  can  make  a  change  of 
variable  A' (  n  )  =  A  ( n )  z  (n)  which  will  produce  an  equation  of  the 

desired  type .  Let  v(n)  be  a  nontrivial  solution  and  the  sum  of 
the  normalizing  series 

S  =  Y.  c  (  k  )  w  (  k )  ^  () 

k=0 


o 


is  known.  Let  N  be  a  large  integer  and  define  zN(n),  0  <  n  < 
N+l  ,  by 

rO  n  =  N+l 

zN(n)  =  < 

ll  n  =  N 

zN(n)  +  a(n ) zN ( n+l )  +  b(n)zN(n+2)  =0,  n  =  N-l.  ....  1.  0 

One  can  approximate  w(n)  by  wN(n) 

wN(n)  =  SzN(n)/SN 

where 

SN  =  E  c(k)zN(k) . 
k=0 

The  algorithm  is  said  to  converge  if 

wN  (  n  )  —  w  (  n  )  as  N  —  oc  . 

The  function  M(a.b:x)  satisfies  the  recurrence  relation 

(2n  +  b+2  )  (  n+a  )  z  (  n  )  -  ( 2n  +  b+l  )  { (  2a-b)  +  +  ^  (x2 ”  +  b+ 2 ) }z  (  n  +  1  ) 

-  ( 2n+b ) ( n+b+ 1 -a ) z ( n+2 )  =  0  . 

The  minimal  sol ut ion  is 

1  n 

(  b  )  2n 


w  (  n  ) 


M ( a  +  n . 2n  +  b : x  ) 


where 


(c)n  = 


T (n+c ) 
F(c) 


The  normalization  relationship  used  in  our  subroutine  is 


S  =  b-1  =  g  ^rrr-  ( b  - 1 )  k  (b+2k- 1 )  w  ( k)  . 

1 _ A  ^  * 


An  obvious  modification  must  be  made  if  b  =  1.  The  algorithr 
is  not  defined  if  b,  b+l-a.  a  are  negative  integers  or  zero. 
The  function  U(a,b;x)  satisfies  the  relationship 


(n+a) (n+a+1 -b) z (n)  -  (  n  +  1  )  [2  (  n  +  a+1  )  +x- b]  z  (  n-t- 1  ) 


+  (n+1 ) (n+2) z (n+2)  =  0 


The  minimal  solution  is 


,  .  xn (a) n (a+1 -b) n  , 

w  ( n  )  =  - - — - -  l  (a+n  ,  b  :  x 


for  |  arg  x|  <  tt  .  A  normalization  relation  i: 


i  =  T.  w(  /  • 

k=0 


In  the  next  section  we  give  a  listing  of  the  Fortran 


su brout i nos . 


1 


Subroutine  Mi  lie r 


SUBROUTINE  M I LLER (N , ALPHA . BETA . X , S . SS . COEFF) 
INTEGER  N 

REAL*8  ALPHA, BETA, X  ,  SS 
REALMS  S ( 0 : 1 OOO ) 

EXTERNAL  COEFF 

C  USES  THE  J.C.P.  MILLER  ALGORITHM  TO  COMPUTE 
C  S  (  0  :  N  )  . 

C  BEGIN 

INTEGER  NN.K 
REAL~8  T , D . EPS , A , B . C 
REALMS  OLDS (0:1000) 

EPS  =  0.000000001 
C  INITIALIZE  OLDS. 

DO  1  k  =  0,  1000 
OLDS (K)  =  0 

1  CONTINUE 

C  CHOOSE  INITIAL  NN . 

NN  =  N  4-  2 

C  INITIALIZE  K,  S  AND  T. 

2  K  =  NN 
S(K+1 )  =  0 
S(K)  =1 

CALL  COEFF (K . ALPHA , BETA , X , A . B , C) 

T  =  2~C~S(K) 

c  TAKE  A  BACKWARD  RECURRENCE  STEP  AND  UPDATE  IT. 

3  K  =  K  -  1 

CALL  COEFF (K . ALPHA , BETA . X . A . B , C) 

S ( K )  =  A~S(K+1)  +  B«S ( K+2) 

C  CHECK  FOR  OVERFLOW  AND  RESCALE  IF  NECESSARY . 

D=  DABS ( S (K ) ) 

IF  (D  .GT.  1.D30)  THEN 
C  BEG  I N 

C A LL  SCALE ( K , NN . S , T . D ) 

END  IF 

IF  (K  .GT.  0)  THEN 
C  BEG  I N 

T  =  T  +  2~C~S(K) 

GO  TO  3 
END  IF 

I  —  1  +  C  *  S  (  0  ) 

DO  4  K  =  0.  N 

S ( K )  =  S(K)/T 

4  CONTINUE 

C  TEMPORARA'  PRINT  STATEMENT. 

C  PRINT*.  S(0) 

C  TEST  FOR  CONVERGENCE . 

D  =  0 

DO  5  K  =  0.  N 
3  D  =  I)  +  S  (  K  )  * » 2 

CONTINUE 
D  =  DSQRT(D) 

T  =  O 


DO  6  K  =  0 ,  N 

T  =  T  +  (S(K)  -  OLDS (K) ) **2 

6  CONTINUE 

T  =  DSQRT(T) 

C  TAKE  ANOTHER  STEP  IF  NO  CONVERGENCE. 

IF  (T  .GT.  EPS~D)  THEN 

C  BEGIN 

NN  =  2*NN 
DO  7  K  =  0,  N 

OLDS (K)  =  S (K) 

7  CONTINUE 

I F ( NN  .LE.  1000)  GO  TO  2 
PRINT  999. NN, ALPHA. BETA, X.T 

999  FORMAT  ( '  **  NO  CONVERGENCE  *«  NN  AP  CP  X  T  ' .  15 . 4E1  *4 . 7  ) 

END  I  F 
SS=S(0) 

RETURN 

END 


SUBROUT I NE  COEFF ( N , ALPHA , BETA . X . A , B . C) 

INTEGER  N 

REALMS  ALPHA. BETA. X. A, B.C 

C  COMPUTES  COEFFICIENTS  USED  BY  J.C.P.  MILLER  ALGORITHM  FOR 

C  A  CONFLUENT  HYPERGEOMETRIC  FUNCTION  M(a,b:x) 

c  SEE  JET  WIMP.  COMPUTATION  WITH  RECURRENCE  RELATIONS 
C  PITMAN  1984  PP .  61-62 

C  BEG  I N 

INTEGER  M.K 
REAL-8 . T . U . V , W 
S  =  2-ALPHA  -  BETA 
T  =  N  +  ALPHA 
M  =  2»N 
U  =  M  +  BETA 
V  =  U  +  1 
W  =  V  +  1 

A  =  (S/W  +  U/X) -V/T 
B  =  (N  +  BETA  -  ALPHA  +  1)-U/T/W 
T  =  1 

IF  (N  .GT.  0)  THEN 
C  BEG  I N 

S  =  BETA  -  1 

DO  1  K  =  1  .  N- 1 

T  =  -T-(l+S/K) 

1  CONTINUE 

T  =  -T»( 1+S/M) 

END  II 
C  =  T 
RETURN 
END 


SUBROUTINE  SCALE (K.N .  S.T.D) 
INTEGER  N.K 
REAL-8  T.l) 

REAL-8  S (0:1000) 

C  BEGIN 

INTEGER  .1 

T  =  T/D 

DO  1  J  =  K .  N 

S  (  J  )  =  S(  J  )/I) 

1  CON  I  1  NT  I  . 

RETURN 


END 


SUBROUTINE  COEFU (N , ALPHA , BETA , X , A . B , C) 

INTEGER  N 

REAL* 8  ALPHA.  BETA.X.A.B.C 

C  COMPUTES  COEFFICIENTS  USED  BY  J.C.P.  MILLER  ALGORITHM  FOR 

C  A  CONFLUENT  HYPERGEOMETRIC  FUNCTION  U(a.b;x) 

C  SEE  JET  WIMP.  COMPUTATION  WITH  RECURRENCE  RELATIONS. 

C  PITMAN  1984  PP .  63-64 

C  BEGIN 

INTEGER  M.K 

REAL*8  S.T.U.V.W 

S  =  ALPHA  +  Q FLOAT (N) 

T  =  S  +  1 .DO 
U  =  S* (T  -  BETA) 

V  =  Q FLOAT ( N  +  1) 

W  =  V  +  1  .  DO 

A  =  ( 2*T  +  X  -  BETA)  *V'/U 
B  =  -  V*W/U 
C  =  1 
RETURN 

END 


Remark:  The  program  that  calls  Miller  must  supply  as  a  last 

parameter  either  COEFF  (for  M)  or  COEFU  (for  U) . 


The  subroutines  are  available  on  a  diskette  from  either  author 


upon  request.  These  subroutines  were  tested  extensively  for 
various  values  of  a,  b  and  x. 

Remark  :  If  the  parameter  is  a  negative  integer,  the  solution 

of  the  differential  equation  is 

y  =  Al,n  (  x )  +  B {  1  n  |  x  |  Ln  ( x )  +  £  dm*™} 

m=0 

where  n  =  - a . 

l.n  ( x )  are  Laguerre  polynomials  whose  coefficients  a| 
sat  i s  f y 


i  -  n  - 


i-i 


i  =  2. 


n 


n 


1 


-  n 


The  coefficients  ,im  satisfy 


■*m+l 


(  m  -  n  )  d|V.  + 


2  ( m  -  n  ) 
m-f  1 


(  m+  1  )  2 


m  =  1 


n-  1 


m 


1 

<11—1  )  2 


°  n 


in 


n 


—  n>-  n  -  1  j 
~  2  1  m-1 

(II 


in  = 


• '  m 


n  -f  1  .  n  +  'J 
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